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Abstract — A blind target detector using the time reversal 
transmission is proposed in the presence of channel correlation. 
We calculate the exact moments of the test statistics involved. The 
derived moments are used to construct an accurate approximative 
Likelihood Ratio Test (LRT) based on multivariate Edgeworth 
expansion. Performance gain over an existing detector is observed 
in scenarios with channel correlation and relatively strong target 
signal. 

Index Terms — Complex double Gaussian; time reversal; de- 
tection; channel correlation; multivariate Edgeworth expansion. 



Despite the practical needs to understand TR detection in 
the presence of channel correlation, the results in this direction 
are scarce. To address this challenge, we propose a blind TR 
detector that admits a general correlation structure between 
the TR channels. A closed-form approximation to the corre- 
sponding likelihood ratio is proposed using the multivariate 
Edgeworth expansion. The approximation is constructed via 
the derived exact moments of the underlying statistics. Numer- 
ical simulations show that the proposed detector outperforms 
the detector in [4] by exploiting the TR channel correlations. 
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I. Introduction 

TIME Reversal (TR) is a waveform transmission method 
that focuses the transmitted energy in dispersive 
medium - the channel [1]. It utilizes channel reciprocity and 
obtains the channel state information by sending a probing 
signal. The backscattered signal is then time-reversed and 
retransmitted. The TR signal is shown to be optimal in the 
sense that the transmission realizes a matched filter to the 
propagation transfer function [1]. The concept of TR was 
originally developed for optical and acoustic applications, 
and it is recently introduced as a detection technique in the 
electromagnetic domain [2^1], where the target to be detected 
is embedded in stationary random multipath scatterers. 

In [2,3], the authors assumed that the multipath channel or 
the channel response signal can be ideally estimated using 
probing snapshots. However, the assumption of a perfectly 
known channel or a noise-free signal may not be realistic 
for practical systems due to e.g. measurement noise. Esti- 
mation accuracy depends on the number of snapshots, which 
is limited by the coherence time/frequency of the channel, 
and the sampling rate of the system [3]. To avoid channel 
estimation, the authors in [4] considered a blind TR detector, 
which utilizes only the distribution of the multipath channels. 
The likelihood ratio test for the TR detector was derived 
assuming statistical independence between the two consecutive 
transmissions. However, this assumption is not valid if the 
transmissions are within the coherence time of the multipath 
channel. Using existing detectors in such a scenario will induce 
performance loss. 
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II. Blind Time Reversal Detection 

We consider blind detection of a point target in the pres- 
ence of multipath scatterers as studied in [4]. The detection 
system sends Q probing signals in the spectral domain at 
the frequencies oj q , q 6 [1, Q\. The sampling frequencies 
are chosen such that each frequency bin is separated by 
the coherence bandwidth of the channel and the spectral 
samples are statistically independent. The multipath channel 
at Lo q induced by the scatterers is modeled by a wide sense 
stationary process. We denote the channels experienced by 
the probing signal and the retransmission as C p (ujq) and 
C r (u> q ), respectively. The channel response of the point target 
is captured by a deterministic response T and the probing 
signal at oj q is denoted as S(u! q ). Note that here we consider 
a general correlation structure between C p (oj q ) and C r (ui q ) 
instead of statistical independence assumed in [4]. As a result, 
knowledge of channel coherence time is no longer required. 
In such a scenario, the detector of [4] suffers performance loss 
as will be shown in Section IV. 

After transmitting the probing signal S(uj q ), we write the 
frequency response as 

Z(0J q ) =(T+ C v (u> q ))S{u> q ) + V p {iOq), 

where V p (ujq) is the measurement noise which is distributed as 
a zero-mean complex Gaussian random variable with Power 
Spectral Density (PSD) cr^. Hereafter, we denote V p (ujq) ~ 
£A/"(0, al) [5]. In this paper, we use a white probing sig- 
nal such that S(u) q ) — \jE s jQ with a transmit power 
E s . The received signal Z(ui q ) is then time-reversed or, 
equivalently, phase-conjugated in the frequency domain and 
scaled to obtain the TR signal, STui^q) — kZ(oj q )*, where 

k 



^s/EjLi l-^^g)! 2 i s a energy normalization factor. 
The value of k is shown to be approximately deterministic with 
a relatively small variance compared to its expected value [3]. 
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The TR signal STn(oo q ) is subsequently transmitted and the 
channel response of the retransmission is calculated by 



Ztr^q) = (T + C r {u> q ))S T R{w q ) + V r (uj q ) 
= X(cj q )Y(u Jq y +V r (uj q ), 



(1) 



where X(u q ) = T + C r (u q ), Y(uj q ) = SrR(u q )* and 
V r {uj q ) ~ CA/"(0,er 2 ) is the measurement noise of the re- 
transmission. In blind TR detection, the channels C p (u! q ) and 
C r (uj q ) will not be estimated by the detector and are only 
known by their statistical distributions. Therefore, a hypothesis 
test can be formulated as follows: in the null hypothesis Ho, 
the target is not present and T = 0; in the alternative 
hypothesis Hi, \T\ > 0. 

We assume the channels C p (uj q ) and C r (ui q ) admit a bivari- 
ate zero-mean complex Gaussian distribution with a common 
PSD P c (uj q ). The correlation coefficient p c between C p (ui q ) 
and C r (uj q ) is defined as p c = E[C p (uj q )C r (uj q )*]/P c (co q ), 
where the notation (•)* is the complex conjugate. In practical 
systems, the channel statistics can be estimated by taking 
snapshots of channel samples and replacing the statistical ex- 
pectation by the sample mean. The measurement noise V p (ui q ) 
and V r are independent of each other and the multipath 
channels. If we ignore the noise term V r (uj q ) in (1), Zt_r(w 9 ) 
is distributed as the product of two complex Gausian random 
variables with 



X(w q )~CM(T,P c (w q )), 




Pc{Wq) 



+ < 



(2) 
(3) 



By definition, the random variables X(co q ) and Y(uj q ) are 
jointly complex Gaussian distributed with a correlation co- 
efficient calculated as 

p= — Pc (4) 

y/l + alQUP^Es) 

To clarify the considered problem, we introduce the random 
variable V(u) q ) = X(ui q )Y(uj q )* and denote its corresponding 
PDF in the complex plane as fr g (Pi,P2;T). The LRT of the 
blind TR detection is calculated by 

i = n 4*^3 £ *o, (5) 



£i fr g (Pi,P2]0) Ho 

lo being a threshold. In the next section, we first derive 
the characteristic function of the product V(co q ). Based on 
this, an exact expression is obtained for fv q (Pi,P2] 0) and an 
asymptotic approximation is given for f-p q (pi 7 p 2 ;T). 

III. Correlated Time Reversal Channel 
A. Characteristic Function 

We first derive the characteristic function ip-p(t) (t 6 
C) of the product V{uj q ) = X(uj q )Y(u> q )* . For notational 
simplicity, the frequency variable uj q is hereafter dropped. 

Recalling equations (2)-(4), the joint PDF of X and Y is 

CXP { fl( ,,, )/( i |,p> 

9{*,V) = ™ + ™ 2 -2?i[p™^] and 



given by [5] as fx, Y (x,y) 
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Pi, Oi (i G {X, Y}) refer to the mean and variance of the 
corresponding random variable. Given Y, X is conditionally 
complex Gaussian distributed with mean px\y = Px + p(y — 
Hy)&x/&y and variance cr^-, y = <7^(1 — \p\ 2 )- Denote 
the real and imaginary parts of V by V\ and V 2 . It is 
straightforward to show that V\ and V 2 conditioned on Y are 
conditionally independent. They follow conditional Gaussian 
distributions 



Vx\Y 
V 2 \Y 



^(5R [ V *px\y]^x\y\v\ 2 /^), 
^(q [ V *px\y] ,a 2 xlY \y\ 2 /2) . 



Therefore, the conditional characteristic function of V given 
Y is expressed as [5] 



ifo\ Y {t\y) = E [exp {m [t*V]} \Y = y] 

= exp{m[t*y*p x \Y} - ^<Tx\ Y \y\ 
The marginal PDF of Y is given by 

f Y (y) = l/(na Y )exp{-\y- f i Y \ 2 /a Y }, yeC. (7) 



2 \t\ 2 }. (6) 



The characteristic function of V can be now obtained by direct 
integration of (6) over the marginal PDF (7) as 



V> 7 



(*) = / ipr\Y(t\y)f Y (y)dy. 

JyEC 



(8) 



Substituting (6) and (7) into (8), we obtain 
ip-p(t) = cxp{-|^ y | 2 /CTy}/(7rcrf,) 



x / exp 

yec 
2^[p Y y*'_ 



°xW 



a 2 Y 4(1 - |pP 



+ 



a 



Y 



4 

.* ( °x 



l^[t*p\ 



ox 



<7Y 

dy. 



(9) 



Applying [6, eq. (3.323/2)] and integrating (9) over real and 
imaginary parts of y, we get 

„ |2^2 _l |,„ ,\2„2 



1 exp/ \^\ 2<J Y + \^\ 2(T x ltl 2 
(*) eP l 4G(t) 



G(t) 

OxOY^x^YP] | t |2 + i&lMxVYt] 



2G{t) 



Git) 



(10) 



where G(t) = 1 + \a\o 2 Y (\ - |p| 2 )|i| 2 - ia x a Y ^[t* p}. 



B. Joint PDF 

Based on (10), we now calculate the joint PDF fv(Pi,P2', 0) 
under the null hypothesis % . When T = 0, p x = Py = 
and (10) is reduced to i>-p(t) = 1/G(t). Applying the inverse 
transform of characteristic function, the joint PDF becomes 

exp{-i&[t*p]} 



fv(Pi,P2;0) 



(2tt) 2 J teC G(t) 



■dt 



2 f 25R[p*p] 
exp <^ 



c 



(11) 



where p = Pi + ip 2 and c = a x & Y (l — \p\ 2 )- Here, the 
function K (-) is the modified Bessel function of the second 
kind [6, eq. (8.432/6)]. The second equality of (11) is obtained 
by using [6, eq. (3.354/5)] and the definition of K (-). 
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Next, we derive an asymptotic approximation to the joint 
PDF f-p(pi,p 2 ;T) using the multivariate Edgeworth expan- 
sion. The Edgeworth expansion was considered as an extension 
to the central limit theorem and developed in the form of a 
moment series expansion weighted by the Gaussian PDFs [7]. 
This method is especially useful when the random variable 
of interest is approximately Gaussian and its moments are 
easy to obtain. For the problem at hand, we can prove 
that fp(pi,p 2 ;T) converges to a Gaussian PDF as fix and fly 
go to infinity. Motivated by this fact, we give the Edgeworth 
approximation to fp(pi,p 2 ;T) based on the closed-form 
expressions for the joint moments of V\ and Vi- 

Let W = (V — fip)/<jp, where \ip = fixPy + P a x^Y 
and cr 2 , = \fix\ 2 a Y + | Py | 2 a\ + a \ a Y- The characteristic 
function of W reads 



^w{t) = exp {-i^R [fipt*] jap) tp v (t/a-p) 



(12) 



We denote 5x = fix I ax and Sy = fiy/ay and notice that 
as \5x\ and |<5y| go to infinity (12) reduces to cxp{ — |t| 2 /4}, 
which is the characteristic function of the standard complex 
Gaussian random variable. Thus, the considered variate V is 
approximately complex Gaussian with mean fip and variance 

Before constructing the multivariate Edgeworth expansion 
for fp(pi,p2]T), we need the following two propositions: 

fix + Vx and Y = 



Proposition 1. Let X 

fiy + Vy, where Vx 

CAf(0,a Y ) with the correlation E[VxV Y \ = pa x ay. 

Denote A tii = {V x , ■ ■ ■ , V x , Vy, ■ ■ ■ , Vy} and B t j = 



CAf(0,a 2 x ) and Vy 



{V x ,--- ,v x ,v, 



Y , ■ ■ ■ 



Vy}. The joint moment M m ,n 



E[P m (V n )*} is given by 



Mr 



m-\-n t t 

EEE 

t=0 i=0 j=0 



n 
t-j 



n 
t - i 



k=l 

where denotes the i-th element of the corresponding 

vector, 7r defines a permutation of the integers 1, . . . , t and fl 
is the set of all t\ distinct permutations. 

The proof of Proposition 1 is a direct application of the 
moment theorem in [8], which is omitted here. 

Proposition 2. The joint moments of V\ and Vi are 

( E[7>rP 2 ] 

E p,m-lpl] 

where J„ +1, the element of the (k+l)-th row and (l+l)- 
th column of matrix 3 m , is given by 





( 


Mm,0 


\ 






Mm-1,1 


- T 1 








/ 


\ 


M QlTn 


I 



LV2J 

j(fe+i, _ ^ i i ~ 2h 

h=0 



m, — 2k 
I - 2h 



(13) 



where [a\ denotes the largest integer less than a. 
Proof: First consider the case m > 2k, 

M m - k . k = ^[{V 1 + iV 2 ) m - 



-k,k 

-2k k 



EE^ 



2k {vl + vl) 

■pm-j 



E 



2h T) j+2h 
' 2 



Summing over I with I = j + 2h, the coefficient of the l- 
th term achieves (13). In the same manner, (13) can be also 
proved when m <2k. ■ 
Following the procedures in [7], the Edgeworth expan- 
sion f s {pi,P2) of the joint PDF fp{p\ 1 p2\T) can be rep- 
resented as 



fp(p!,p 2 ;T) 



<P{pi,P2) + ^L 



{-4>;x v ){puP2), (14) 



where (j)(pi,p 2 ) is the bivariate Gaussian PDF with mean 
H = [fip 11 fip 2 ] T and co variance matrix R of V\ and TV 
Here x„ with v = \v\,v 2 \ refers to the joint cumulant of 
V\ and Vi- These can be readily calculated by the mapping 
between joint cumulants and moments, which are computed 
by Proposition 1 and 2. It is noted that the expansion (14) 
uses up-to s-th order joint cumulants of V\ and V2 and 
matches the first (s — 2) moments between V and the ap- 
proximation (14). The term Lj (—<p;xv) (PI1P2) defines the 
Cramer Edgeworth polynomial Lj(ti,t 2 ; Xv) w i tn ea ch t^t^ 2 
replaced by H(pi,p 2 ;v,'R~ 1 )(j)(pi,p 2 ), where H denotes the 
multivariate Hermite polynomials [9]. The Cramer Edgeworth 
polynomials Lj{ti,t 2 \Xv) can be generated by the formal 
identity between two power series [7] 



^2L j (t 1 ,t 2 ;Xv)u r = 

r=l 

where x«(*i,t2) = E„ 1+V2 = s l^k 




XvtYt. 



2 ■ 



IV. Numerical Results 

In this section we investigate the accuracy of the derived 
approximation to fp(pi,p 2 ;T) by simulations. Then we com- 
pare the performance of the proposed detector with the one 
studied in [4]. 

A. Approximation to fp(pi,p 2 ',T) 

In Fig. 1 we plot the approximative PDF (14) and compare 
with the simulated fp(pi,p 2 ;T) assuming fix = 2 + 2.5i, 
fiy = 2.1 + 1.8i, a x = <y Y = 1 and p = 0.3 + 0.3i. 
Fig. 1 shows that the Edgeworth approximation achieves a 
good agreement with the simulation using up-to 6-th order 
joint cumulants of V\ and V 2 . The approximation is less 
accurate near the origin since the function fp(pi,p 2 ;T) has 
a simple singularity at the point (0, 0). 

To further illustrate the accuracy of the approximation, in 
Table I we tabulate the Mean Square Error (MSE) of the 
approximative PDF (14), defined as r {a , P2}eK 2_|/ s (pi,P2) - 

fv{Pi,P2', T)\ 2 &F v (pi,p 2 ; T), where fp and Fp are empir- 
ical PDF and CDF of V. In addition to the MSE calculated 
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Fig. 1. Plot of two-dimensional PDF f-p(p\,p2\T) for px = 2 + 2.5i, 
py = 2.1 + 1.8i, ax = °V = 1 and p = 0.3 + 0.3i. (a) approximation (14) 
with s = 6 and (b) simulation. The mean square error of approximation is 
9.96 X 1(T 7 . 

TABLE I 

Mean square error of approximative distribution (14) 





1 + 1.25i 


2 + 2.5i 


4 + 5t 


p Y 


1.05 + 0.9i 


2.1 + 1.8* 


4.2 + 3.6i 


MSE 


2.78 x 10" 4 


9.96 x 10 _v 


2.99 x 10~ 8 



with the parameters used for Fig. 1, we also calculate the 
MSEs with both nx and fj,y having 3 dB decrease and 
increase, respectively. The results show that the accuracy of the 
approximation (14) is improved as the magnitudes of fix and 
Hy increase, which is in line with the analysis in Section III-B. 
This conclusion has been also verified with other parameter 
combinations. 

B. Performance Comparisons 

Here we consider a scenario where the point target has 
a constant response T = e l7r / 4 and the multipath chan- 
nels are of equal PSD P c (w 9 ) = P c over the frequency 
bands {u q }. The signal-to-clutter ratio and signal-to-noise 
ratio are defined as SCR = 101og 10 (|T| 2 / ' P c ) and SNR = 
101og 10 (E S \T\ 2 /{Q<tD), respectively where E s = 1. The 
performance of the proposed LRT detector (5) under channel 
correlation, denoted as LRT-C, are evaluated with correlation 
coefficients p c = 0.1 + OAi and 0.1 + 0.71 In addition, we 
consider the cases where the target has a relatively strong 
channel response (SCR = 5 dB, SNR = 5 dB) with sample 
size Q = 5, denoted by the square markers in Fig. 2 and 
relatively weak target response (SCR = dB, SNR = dB) 
with Q = 20, denoted by the circle markers. 

Fig. 2 shows results from Monte Carlo simulations for the 
Receiver Operating Characteristic (ROC), where the detection 
probability (Pd) is plotted as a function of the false alarm 
probability (Pfa). The ROC of the proposed LRT-C detector 
is evaluated by test (5). For comparisons, we also computed 
the ROCs of the LRT detector designed for independent TR 
channels (LRT-I) [4]. From Fig. 2 we can observe that the 
proposed LRT-C detector outperforms the existing LRT-I by a 
substantial margin when the target is relatively strong. These 
results are expectable since the LRT-C detector utilizes the 
channel correlation, which increases the extent of coherence 
between the TR signal and the multipath channel. As p c in- 
creases to 0.1+0.7i, the blind TR detection tends to a coherent 



(a) p c = 0.1 + OAi 




Fig. 2. Receiver operating characteristic for the proposed LRT-C (5) and 
LRT-I in [4] with relatively strong target signal (SCR = 5 dB, SNR = 5 dB) 
with sample size Q = 5, denoted by and relatively weak target signal 
(SCR = dB, SNR = dB) with sample size Q = 20, denoted by 'O'- 
(a) p c = 0.1 + OAi and (b) p c = 0.1 + 0.71 

detection and the proposed LRT-C detector achieves improved 
detection probability at a fixed false alarm probability by 
capturing the channel correlation. When the target response 
is weak with small channel correlation, the performance of 
LRT-C becomes worse than the LRT-I. This observation is 
consistent with the Neyman-Pearson theorem [10], as the 
increased approximation error of (14) under this condition 
leads to a non-trivial deviation from this optimal detector. 

V. Conclusion 

We proposed a blind time reversal detector which works 
in the presence of correlated channels. Using Edgeworth 
expansion, a simple and accurate closed-form approximation 
was derived for the likelihood ratio test. Simulations show 
the superiority of the proposed LRT-C detector in scenarios 
with arbitrary channel correlation and relatively strong target 
response signal. 
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